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Abstract 

Many calculations for the production of light and intermediate particles resulting from heavy ion 
collisions at intermediate energies exist. Calculations of properties of the largest fragment resulting 
from multifragmentation are rare. In this paper we compute these properties and compare them 
with the data for the case of gold on carbon. We use the canonical thermodynamic model. The 
model also gives a bimodal distribution for the largest fragment in a narrow energy range. 
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I. INTRODUCTION 



The statistical model of nuclear disassembly in heavy ion collisions is quite successful. 
Here one assumes that the disintegrating system with a given excitation energy expands 
to greater than normal volume before it breaks up into many composites of various sizes. 
Interactions between different composites in the rarefied situation can be neglected and 
the break up can be calculated using the laws of equilibrium thermodynamics. This basic 
assumption is implemented in different versions according to degrees of sophistication and 
detail. Thus we have the statistical multifragmentation model (SMM) of Copenhagen the 
microcanonical model of Gross [2] and Randrup and Koonin |3]. An easily implementable 
canonical ensemble model was later introduced [4|. The details of the model and many 
applications can be found in a recent publication 5|. 

It has been customary to pay a great deal of attention to the production cross-sections of 
intermediate mass fragments and light charged particles. Here we study the properties of the 
largest fragment that emerges in multifragmentation. These properties are not easy to study 
but the canonical model allows for such computation. As we will see, the heaviest fragment 
also reveals some interesting physics. Data on the heaviest fragment and fragments with 
the maximum charge were published by EOS collaboration 0, Q, Q]- In these experiments 
one studied the disintegrations of projectile- like excited fragments in the reactions Au+C, 
La+C and Kr+C. We do our calculations for Au+C Calculations for the other two systems 
will be similar. 

In section II we write down the theoretical formulae needed for the calculations. Results 
are presented in section III. Summary and conclusions are presented in section IV. 

II. FORMULAE IN THE CANONICAL MODEL 

The formulae used in this calculation are provided here. Simpler formulae for a hypo- 
thetical system of one kind of particles are given in 5|. 

We will consider disassembly of the projectile-like fragment (PLF) where the PLF is 
formed by the shearing off of a part of 197 Au (by the 12 C target). This PLF will have 
a charge Z (usually less than 79) and a neutron number N (usually less than 118). This 
PLF will break up into many composites with charges i and neutron number j (for example, 
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9 Be has i = 4 and j = 5). If the number of composite with proton and neutron numbers 
i and j respectively is n it j then conservation of charge and baryon number dictates that 

Z = ^2 i x n i,j and N = j x n i,j- 

The canonical partition function for the system at a given temperature T is given in our 
model by 

o^ = EII^r w 

Here the sum is over all possible channels of break-up (the number of such channels is enor- 
mous) which satisfies the conservation laws; cjjj is the partition function of one composite 
with proton number i and neutron number j respectively and rijj is the number of this 
composite in the given channel. This is a low-density/ high temperature approximation the 
justification of which is demonstrated in Q and fl. For a given channel the sum . n^j is 
called the multiplicity of the channel. The one-body partition function is a product of 
two parts: one arising from the translational motion of the composite and another from the 
intrinsic partition function of the composite: 
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(2vrm(z + ])Tf' 2 x z hJ {mt) (2) 

Here m(i + j) is the mass of the composite and Vf is the volume available for translational 
motion; Vf will be less than V, the volume to which the system has expanded at break up. 
We use Vf = V — Vq , where Vq is the normal volume of the PLF. We will shortly discuss 
the choice of Zij(int) used in this work. 

The probability of a given channel P(njj) = P^nQ^,^^,^^ ) is given by 

vz,v Tkj } - 

The average number of composites with i protons and j neutrons is seen easily from the 
above equation to be 

(thj) = "is—?] ( 4 ) 

WZ,N 

The constraints Z = ^ i x mj and N = ^ j x n^j can be used to obtain different looking 
but equivalent recursion relations for partition functions. For example 

Qz,N = 7F i&i,jQz-i,N-j (5) 



Instead of labelling partition functions by their charge and neutron numbers Z and iV we 
can for example label them by total mass (A = Z + N) and charge Z. Labelling them by Q 
we have 

Qa,Z = -J ^ atia,iQA-a,Z~i (6) 

These recursion relations allow one to calculate Qz,n (equivalently Qa,z) very quickly in the 
computer. 

We list now the properties of the composites used in this work. The proton and 
the neutron are fundamental building blocks thus zi^{int) = zo^(int) = 2 where 2 
takes care of the spin degeneracy. For deuteron, triton, 3 He and 4 He we use Zij(int) = 
(2sij + 1) exp(— (3eij(gr)) where (3 = 1/T, eij(gr) is the ground state energy of the compos- 
ite and (2sij + 1) is the experimental spin degeneracy of the ground state. Excited states 
for these very low mass nuclei are not included. For mass number A = 5 and greater we use 
the liquid-drop formula. For nuclei in isolation, this reads (a = i + j) 

z^iint) = exp hw a - a{T)a 2 ^ - k-^ ~ s^—^- + — ] (7) 

The derivation of this equation is given in several places |l, 5] so we will not repeat the 
arguments here. The expression includes the volume energy, the temperature dependent 
surface energy, the Coulomb energy, the symmetry energy and contribution from excited 
states since the composites are at a non-zero temperature. For (the proton and the 
neutron number) we include a ridge along the line of stability. We have used two sets to test 
the sensitivity of the results to the width of of the ridge. We call these set 1 and set 2. In set 
1 for each a between 5 and 40, 5 isotopes are included (case for a=4 and lower was already 
mentioned before); For a > 40, 7 isotopes are included for each a. In set 2, 5 isotopes are 
used for a between 5 and 9, 7 isotopes are used for a between 10 and 40 and 9 isotopes for 
each a > 40. The results are quite similar in most cases except for one figure which we will 
point out when we compare with data. It should be pointed out that enlarging the width 
of the ridge does not necessarily imply a better calculation as one may begin to overcount 
the phase space. 

The Coulomb interaction is long range. The Coulomb interaction between different com- 



posites can be included in an approximation ca 
incorporate this following the scheme set up in 



led the Wigner-Seitz approximation. We 
l|. This requires adding in the argument 



of the exponential of Eq.(7) a term k^/w (K)/^) 1/3 (^^) 1/3 . Here A , Z are the mass and 
charge number respectively of the disintegrating system, V is the normal nuclear volume 
for this system and V, the freeze-out volume (typically 4-5 times Vo). Defining i? 3 = 

3.Z 2 

the energy of the system is given by = + z2i j( n i,j) e i,j where for a > 4 we have e^j = 
%T+a(-W +T*/eo)+a(T)a 2 ' 3 +K^[1.0-W 

For a < 4 we use e id = §T + e^^r) - K^(Vb/F) 1/3 (^^) 1/3 . We label as E* the excita- 
tion energy: E* = E — E(gr) where E(gr) is calculated for mass number A and charge Z 
using the liquid-drop formula. 

The central issue in this work is the calculation of (Z max ) and (A max ) (and their fluc- 
tuations) when the fragmenting system has charge Z and mass number A . Let us fix 
on calculating (Z max ). The calculation for (A max ) can be done by analogy with appro- 
priate and obvious changes in the formulae. There is an enormous number of channels 
in Eq.(l). Different channels will have different values of Z max . For example there is a 
term -§f-^f in the sum of Eq.(l). In this channel Z max is 1. The probability of this 
channel occurring is (from Eq.(3)) qJ~^~zT ~nT- The full partition function can be writ- 
ten as Qz,n = Qz,n(^i,o, Wo,i, ^1,1) ^ij )• If we construct a Qz,n where we set 

U! Z U) N 

all u;'s except u^o and ouo,i to be zero then this Qz,n(wi,o, wq,i, 0, 0, ) = ^fivf 

and this has Z max = 1. Consider now constructing a Q(Z,N) with only three cj's: 

<5z,at (<^i,o, ^0,1,^2,1, 0, 0, 0, ). This will have Z max sometimes 1 (as -§f-^f is still there) 

and sometimes 2 (as, for example, in the term -ft^ prj^i ^ )■ 

We are now ready to write down a general formula. Let us ask the question: what is 
the probability that a given value z rn occurs as the maximum charge? To obtain this we 
construct a Qz,n where we set all values of u z j = when z > z rn . Call this Qz,N(z m )- 
Then Qz,N(z m )/Qz,N (where Qz,n is the full partition function with all the cj's) is the 
probability that the maximum charge is any value between 1 and z m . Similarly we construct 
a Qz,N(z m — 1) where u z j is set at zero whenever z > z m — 1. The probability p that z m is 
is given by 



n (y \= n (7 —7 \ 9^A^A ~ Qz,N(z m - 1) 




(8) 



The average value of Z max at given temperature and for given Z , A is 




(9) 
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and the fluctuation is 

RMS(Z max ) = \J^2p(z m )(z m - Z max ) 2 (10) 

Before we end this section we want to mention that the largest composite we obtain at 

the end of the above calculation is at a finite temperature T and can further decay by 

evaporation, ending up in a lower mass or charge number. From our past experience we 

know what the effect will be: Z max will decrease slightly but the effect on the RM S (Z max ) 

will be significant. Without the evaporation the calculated RMS will be overestimated (see 
fl 

Fig. 19 of [5]). Inclusion of evaporation as an after burner is not possible at this stage. The 
sum over z m (Eq.(9)) is too huge and for each z m there is a sum over neutron number. We 
will hope to get Z max nearly right but will overestimate the RMS value. The next section 
compares data with our calculations. 



III. DECAY OF EXCITED PROJECTILE-LIKE FRAGMANT 

In the EOS experiment part of the projectile (Au, La and Kr) is sheared off by the 12 C 
target. We will compare our calculations with the Au+C data. The size of the excited 
PLF which decays depends upon the impact parameter which also determines the amount 
of excitation energy per nucleon in the excited PLF. In Figs. 1(b) and 1(c) of [6| the size of 
the excited PLF (the size and charge are denoted by A and Zq) is plotted as a function of 
excitation energy E* per nucleon. This aspect of the experiment depends upon dynamics 
and is outside the scope of a thermodynamic model. However, given E* and Aq, Zq, this 
excited PLF will expand and break up into many pieces and this is calculable in a canonical 
thermodynamic model. In our Fig.l (data taken from Fig. 1(f) of [6]) we have plotted 
Zmax/Zo as a function of E* per nucleon where Z max is the average maximum charge carried 
by a composite. The data are shown as points joined by a dotted curve. The other curves 
are exploratory calculations. Our fit to experimental data is given in Fig.2. We now explain 
how the calculations are done. 

For a given value of Aq, Zq, the experiment provides the value of E* . The beginning of a 
canonical thermodynanic model besides A , Zq, are a temperature and a freeze-out density 
p/po (freeze-out density in unit of normal nuclear density po= 0.16fm -3 ) which will then 
provide all the observables including E*. For a fixed E* the canonical model employs a 
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fixed p/po- For example, for central collisions of Sn on Sn at 50 Mev per nucleon beam 
energy, a freeze out density of one-sixth normal density gives good results for production 
cross-sections of intermediate mass fragments ^]. However, in the EOS experiment E* 
varies over a wide range (from less than 2 MeV per nucleon where the validity of the 
thermodynamic model as used here can be questioned to 10 Mev per nucleon where the 
thermodynamic model is expected to work well) and hence we should expect that p/po 
will also need to vary in this interval. In general, the freeze-out density will decrease as 
M/Aq (M=multiplicity, A =mass number of the dissociating system) increases, reaching 
some asymptotic value for large multiplicity. In SMM [l| the freeze-out density varies in 
each channel, decreasing as the multiplicity increases (this makes Monte-Carlo simulation 
mandatory). In the canonical model, at a given temperature, the freeze-out density is kept 
fixed irrespective of channels. Thus the freeze-out density can be dependent only on the 
average multiplicity. Past comparisons with SMM predictions showed that at least for the 
obsrvables studied so far this simplification is quite adequate [lol ]. 

Fig.l compares the data if in the calculation the freeze-out volume is kept fixed at p/po 
=0.25 (a typical canonical model value) or at 0.39 (this is an often quoted value in the 
model used in The value 0.39 is clearly better at low values of E* but leaves too 

large a residue at higher E* whereas the value 0.25 is better at the higher end of E* but 
is an underestimation at the lower end of E*. The data undoubtedly point to the need 
of a variable p/po if the whole spectrum of E* is to be covered. In Fig.l, for brevity we 
show results with set 1 only (see previous section for the range of nuclei covered in set 1). 
Fig. 2 compares data with our calculation where we use a variable p/po- We have used a 
parametrisation p/po = a + 6exp(— c(E*/A )) where a=0.17, b = 0.83 and c = 0.417 MeV -1 . 
No optimisation of the fit was tried but the values are suitable for the low and high limits 
of E*. Our values of p/po are also very similar to those quoted in Table II of jsj where a 
different model was used. We have shown results with set 1 and set 2 (this includes a larger 
number of composites than set 1; see previous section). 

In the canonical model calculation normally the inputs are the freeze-out density and 
temperature. Here we use the freeze-out density and E*. Given this density and E* we find 
the temperature which would give back the E* we started with. We then calculate Z max . 
Calculations for A max were also done and the fits are similar. 

Next we turn to results for 72 and RMS(Z max /Z Q ) which are shown in Fig. 3. Data are 
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from Fig.2 of [6]. Here 72 = M j^i° where is the fcth moment of the fragment distribution: 
Mfc = ^2,n a a k . The over-estimation in the calculation for RMS was already alluded to in 
the previous section. Afterburner for the largest fragment will bring the distribution more 
closely packed near the line of maximum stability thus reducing the fluctuation. 

The probability distribution (Eq.8) of Z max as a function of Z max for EOS experiments 



is not known to us but in view of the recent interest in such distributions [11|, [12] we have 
computed the distributions and displayed them in Fig. 4. We find that there is a window 
where the distribution is bimodal. Usually for discussion of bimodality one uses a fixed 
freeze-out volume and varies the temperature but here, as we discussed before, our freeze-out 
density actually changes as E* (hence T) changes but the nature of bimodality is still clearly 
seen (if we fix the freeze-out density the plots are very similar). Probability distribution for 
Amax (calculated but not shown here) also has the window for bimodality. Connection 
between bimodality and first order phase transition in the canonical model (which does 
have a first order phase transition) is being pursued and we also hope to do calculations for 
other experimental data. 

IV. SUMMARY AND DISCUSSION 

The canonical thermodynamic model clearly reproduces many important features of the 
largest fragment resulting from multifragmentation. The model is very easy to implement 
and yet is very realistic. The same model gives remarkable fits to light charged particles, 
intermediate mass fragments and as we have just seen, the largest fragment as well. We 
have thus far reported the calculation of properties of the largest fragment without making 
any link with the subject of liquid-gas phase transition. We make a brief connection here. 
In a large system(see details in a recent article [131 ]) composites with charge z < 20 comprise 
the gas phase. At co-existence, there will be, in addition, one large composite which is 
the liquid. This picture gets blurred as the system size is reduced, but , nonetheless, the 
largest fragment (at temperatures below that which displays bimodality) is an approximate 
scaled down version of the liquid. In a finite system bimodality in the largest fragment 
distribution occurs in the temperature (energy) window where the system passes from the 
liquid-gas co-existence phase to the pure gas phase. The largest cluster has been discussed 



in the literature before 



fl Q, Q. 



161 ] . However these all use the lattice gas model or the 
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Ising model (often with fixed magnetisation). Our approach here is more directly related to 
nuclear physics phenomenology with phase transition issues in the background. 
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FIG. 1: Experimental points (from Fig. 1(f) of [a] ) are joined by a dotted line to guide the eye. The 
calculations are exploratory: these are done for a fixed freeze-out density (solid line for p/po=0.25 
and dashed line for p/po= 0.39). A higher freeze-out density produces a higher value of (Z max ). In 
the figures Z max stands for this average value. The calculations shown here include the composites 
of set 1. See the discussion following eq.(7) for enumeration of composites included in set 1 and 
set 2. 
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FIG. 2: It is expected that the freeze-out density will decrease with excitation energy in this energy 
range. The calculations here use a parametrisation p/po = a + 6exp(— c(E*/Aq)) where a=0.17, 
6=0.83 and c=0.417 MeV -1 . The solid line uses set 1, the dashed line uses set 2. Experimental 
data are joined by a dotted line. 
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FIG. 3: Top panel compares experimental data for 72 (fig. 2 of [6]) with calculations. For definition 
of 72, see the text. The solid curve uses set 1 and the dashed curve uses set 2. The bottom panel 
compares experimantal values of RMS of Z max with the canonical calculation. As explained in the 
text, without an afterburner evaporation calculation, RMS values from theory will significantly 
overestimate actual RMS. Evaporation will bring the population much closer to the line of maxi- 
mum stability. For an estimate of effects of evaporation see fig. 19 of [5]. It is not surprising that 
with set 2 (dashed curve) calculated values of RMS exceed those calculated with set 1. As more 
nuclei further from the line of maximum stability are included the effects of evaporation will be 
stronger. 
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FIG. 4: Theoretical values of the probablity of occurrence of a Z max plotted as a function of Z max 
(see eq.(8)). In a small window of E* (equivalently T as shown here) a bimodal distribution is 
seen. Similar plot of probability of A max against A max also shows bimodality. The plots here use 
set 1 
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